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Abstract 


This paper presents a numerical solution for a time delay parabolic 
problem (reaction-diffusion) containing a small parameter. The numerical 
method combines the implicit Crank—Nicolson scheme for the time deriva- 
tive on the uniform mesh and the central difference scheme for the spatial 
derivative on the Shishkin-type meshes. It is shown to be second-order 
uniformly convergent in time and space. Then Richardson extrapolation 
technique is applied to enhance the accuracy from second-order to fourth- 
order. The error analysis is carried out, and the method is proved to be 
uniformly convergent. These two methods are applied to two test examples 
in support of the theoretical results. 


AMS subject classifications (2020): 65M06, 65M12. 


Keywords: Time delayed parabolic problem; Boundary layer; Post process- 
ing technique; Singular perturbation. 


* Corresponding author 
Received 7 August 2021; revised 10 November 2021; accepted 20 November 2021 
Jugal Mohapatra 


Department of Mathematics, National Institute of Technology Rourkela, 769008, India, 
e-mail: jugalQnitrkl.ac.in 


Lolugu Govindarao 

Department of Mathematics, Amrita School of Engineering Coimbatore, Amrita 
Vishwa Vidyapeetham, Amrita University, Ettimadai, Tamilnadu, India, e-mail: 
l_ govindarao@cb.amrita.edu 


250 


A fourth-order optimal numerical approximation and its convergence ... 251 


1 Introduction 


In singularly perturbed delay partial differential equations (SPPDEs), a 
small parameter affects the highest derivative of partial differential equations 
(PDEs) and also contains more than one delay term in the time direction. 
These SPPDEs have existence in subjects like physical/chemical science, bi- 
ology, and ecology. Time delay (large) diffusion problems (the present time 
depends upon the past) appear in mathematical models in population dy- 
namics [13, 26, 17] and also in biological modeling [29]. 


Consider the following time delay reaction-diffusion initial-boundary- 
value problems (IBVPs): 


(5 +£e5)2(y. 9) =—by, s)2ly, s—D+fluss), (ys) €D, 


Os 
z(y, 8) = Oo(y, s), (y, s) € 0D, (1) 
2(0,s) = O1(s), on 0, = {(0,s):0<s<S}, 
z(1,s) = 0,(s), on 6, ={(1,s):0<s<S}, 


where Leyz(y, $) = —eZyy(y, $s) + p(y)z(y, s). Here = (0, 1), D =H x 
(0, S], 0 = 0,U 8,U5,. Moreover, 7 > 0 is a given constant. This model (1) is 
singularly perturbed, parabolic in nature with 0 < ¢« < 1. Moreover, 6, and 
d; are the right and the left sides of the domain D, and 6, = (0, 1] x [—7, 0]. 
The functions b(y, s) (b(y, s) > 0), p(y)(p(y) = B > 0), the source term 
f(y, s) on D and O/(s), ©,(s), Ox(y, s) on , are bounded, sufficiently 
smooth functions. The assumption on the time S is S = my for m € N. 
We also assume that suitable compatibility conditions hold true on the given 
boundary and initial data to ensure a unique solution for (1) that exhibits 
boundary layers along with end points of y [1, 19]. These conditions are 
©,(0,0) = 07(0), E,(1,0) = O,(0), and 


d@;(0) 07@,(0, 0) 
ds © Oy? 

d@,(0) _820,(1, 0) 
ds 7 Oy? 


p(0)©2(0, 0) = —b(0, 0)02(0, = 7) + F(0, 0), 


+ p(1)O,(1, 0) = —b6(1, 0)O,(1, — 7) + fC, 0). 


Due to the small parameter ¢ involved in the problem (1), the classical 
methods on equal step length for solving (1) fail to give accurate results. They 
are mostly unstable and unacceptable [22, 27]. Hence, we choose the fitted 
mesh idea through the nonuniform mesh as given in [14, 20, 22, 25, 27] and 
the references therein. One can refer [1, 7, 10, 9, 11, 18, 16, 23] for some order 
enhancing methods of IBVPs. Some articles are available, which discuss both 
the analytical and the numerical techniques for SPPDEs in literature. Ansari, 
Bakr, and Shishkin [1] numerically solved the problem (1) on the Shishkin 
mesh. Das and Natesan [6] gave details of numerical results for the time delay 
parabolic convection diffusion problem. Indeed, the methods discussed above 
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using numerical techniques are of first- or second-order accurate. Hence, order 
enhancing techniques are needed for the problem (1), which is the main aim 
of this work. 

Richardson extrapolation through averaging the numerical solutions com- 
puted on two embedded meshes provides a good approximation to the exact 
solution. This helps to increase the accuracy and the order of convergence. 
Mohapatra and Natesan [21] used the extrapolation technique for solving sin- 
gularly perturbed delay BVPs while Shishkin and Shishkina [28] applied on 
time dependent IBVPs of reaction-diffusion type. This technique is used in 
[5] for convection-diffusion singularly perturbed parabolic problems on the 
adaptive mesh. This article aims to get a fourth-order accurate solution for 
(1) using the extrapolation technique. Initially, the central difference scheme 
is used on Shishkin-type meshes and the Crank—Nicolson method on tempo- 
ral direction on uniform mesh. Here, problem (1) is solved with M number of 
mesh points in spatial and N number of mesh points in time direction. After 
that (1) is solved by the above methods with 2M and 2N number of mesh 
points. Then the rate of convergence increases from second- to fourth-order 
globally by taking a proper combination of these two solutions. 

The rest portion is arranged as follows. In Section 2, we describe the 
continuous solution to the problem. Section 3 studies the construction of the 
numerical schemes. In Section 4, we implement the post-process ideas and 
the theoretical analysis. Section 5 presents the numerical results through 
plots and tables. We denote C and the subscripted C’s as constants, which 
are positive, independent of the small parameter (€) and spatial and time 
mesh sizes. The error is represented in the supremum norm (|| - ||,o). It is 


defined as ||5||.. = sup |6(y,s)| for any function § on the domain D. 
(y,s)ED 


2 Analytic solution and its behavior 


As the perturbed parameter ¢ — 0 in (1), the problem reduces as given below: 


7s + p(y)zo(y, s) = —bly, s)zo(y,s— 7) + fly, 5), (y, 8) € 9D, 


Ozo(y, 8) 
zo(y, s) = O(y, 8), (y, s) € dp. 


(2) 
It is clear that the solution to (1) has boundary layers on 0; and 0,. The 
characteristics of (2) are the vertical lines y = C, which implies that boundary 


) 
layers arising in the solution are parabolic type. The operator (ae + Leg) 


in (1) satisfies the following lemma known as the maximum principle. 


Lemma 1. Suppose that U(y, s) € C°(D)NC?(D) satisfies (+824) Hu, ay > 
S 
0 in D and that U(y, s) > 00nd. Then W(y, s) > 0 for all (y, s) € D. 
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Proof. The proof of this lemma is available in [1]. 


2.1 Solution decomposition 


The continuous solution z(y, s) to (1) is decomposed as z = v, + vs. The 
regular component v, expresses as u(y, $8) = Uro(y, 8)+EUri(y, §), (y, 8) € 
®, where v;9 and v;, are the solutions to the following problem: 


(vro)s(y, 8) + P(y)ro(y, 8) = —b(y, 8)rro(y.s— 7) + f(y, 8), (Y, 8) ED, 
Vro(Y, s) = Os(y, 8), (y, s) € 0p, 


3) 
(5; + Seu)era(y, 8) =—b(y, 8)0ra(y,s—7) + (rola (Ys 8) €D, 


Uri(y, 8) =0, (y, s) ET. 


(3) 
The regular component v,(y, s) satisfies the following problem: 
3) 
(5, + few)ory, 8) =—bYy, svy,s—7) +f, 8), 8) €9, 
Urly, s) = (y, s), (y, s) € Db, (4) 
v,(0,t) = v-9(0, 8), on 0;={(0, s):0<s<S}, 
vr(1,t) = vro(1, 8), on 6,={(1, s):0<s<S}, 
and the singular component v,(y, s) satisfies the PDE 
a) 
(= + Ley) voy, s) = —b(y, s)us(y, So Y) (y, 8) € D, 
Us(y, s) = 0, (y, s) ETs, (5) 
us(0, s) = ©; — vr (0, 8), on 6, ={(0, s):0<s<S}, 
us(1, s) = 0, — vro(1, 8), on 6, ={(1, s):0<5<S}. 


Now, we can write vs, = Usi + Usp, where vs; is the boundary layer part on 0; 
and vs, is the boundary layer part on 0,. Here v,; and vs, are defined by 


(sy t Lew) railys 8) =, s)ouly, 8-7), (8) €9, 


vaiXO, s) = - Uro (0, s) (y, s) € 0, Vet (Y, s) = 0, (y, s) € do U o,, 
az t Lew) Yer(s 8) =—b(y, 8)vsr(y, 8-1), (Ys 8) € 9, 
Usr(1, 8) =O, — Upo(1, 8), (YS) € Oy Usr(y, 8) =9, (y, 8) € WU dp. 
(6) 


Theorem 1. For all integers 1 > 0 and m > 0 with 0 <1+2m < 8, vu, and 
vs, defined in (3) and (5), respectively, satisfy 
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| ie Scie  ). Ge Oeg, 
and 
| ae z < ce“/?( exp(—y B/e) +exp(—(1 — y)V8/2)), (y, 8) €®, 
| aan . < Ce“! ?(exp(—y Be), (y, 8) €D, 
| saat 7 < ce“/?( exp(-(1 — y)V/B/e)), (y, s)€ 9. 


Proof. One may refer [1] for the details. 


3 Discretization 


On [0, S], the uniform time step As is used in the time direction such that 
1 = fao=— nas, A= Ooh. aes, Ae a=oiN, We —4e,. = 
jAs, j =90,...,p, Sp=y, As = y/p}. Here, p represents the number of 
mesh points in [—y, 0] and N denotes the number of mesh points in temporal 
direction on the interval [0, S]. The step size (As) satisfies pAs = y, where 
p> 0 is an integer s, = nAs, n> —p. 


3.1 Discretization of the spatial domain 


1 
Let € = min { me poVeln MW \, where po > 2/2 is a mesh transition parameter. 
We divide II = [0,1] into three subdomains as H = II; U I. UU,, where 
TI; = [0, €], 1. = (€,1 — €], and I, = (1 — €, 1]. Without loss of generality, we 
assume that M is even and that M > 8. 


Shishkin mesh(S-mesh) 


One can find the construction of the S-mesh in [22, 27] briefly described as 
follows: let us define S-mesh as II} = {y; € [0,1],0 <i < M}, where 
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fe 
“ for O0<i<M/4, 
2i(1 —2 M M 
Yi = 2i(1 — 26) for A ee 
M 4 4 
Ai 3M 
“ for -+1<i<M. 


1 
If €= rm" then the mesh has equal step length and otherwise when € = 


poelnM, the mesh is changing at near the end of 0; and 0,, where 
Yi — Yi-1 = 4€M~—'. Therefore, the mesh is piecewise uniform. 


Bakhvalov-Shishkin mesh (B-S-mesh) 


The idea of constructing the B-S mesh is available in [4, 27]. The mesh is 
constructed as follows: 


2 ; 
3\-m (1-20-44) J, for 0<i<M/4-1, 
_ YM/441 —) . M Zee 3M 
Yi = yaya + ( M/2+2 (i — M/4+1) for re ee 
2 1) M-i 3M 
1 2( in (1-201 - 4) 47 )) for -+1<i<M. 


We define the numerical domain 9” = Ng xT on D and 6” = TY x Te 
on 0. 


3.2 Semidiscretization 


The Crank-—Nicolson scheme for the time variable of (1) is given by 
z9 = Op(y, —s;) for j=0,...D, ye, 
T+ Ae Ley gntl — = — prtlizn—ptl _ gnyn-p | fret i) ( Ao.) 2" 
2n+1(0) > O1(sn41), gnrtl(1) = ©,(Sn41), 


(7) 

where f” = f(y,tn), 6” = bly, sn), 2” = z(y, Sn) is the semidiscrete 

approximation to 2(y, s) of (1) at s, =nAs. Let e™t! = 2"t!— Z"*! be the 
local truncation error of (7) and let Z"* be the solution to 


-7 = ©,(y,—t;) for j=0,...p, 2€D, 
1+ A8e.,)artt = $8 (— prtigrop+t — prance 4 pret tym) 4 (7 288.4)", 
2"+1(0) = O1(Sn41), 2" 71(1) = Or(Sn41). 
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Hence, the global error at s” is given by E” = z(y,s”) — z"(y). 


3.3 Bounds on the solution and its derivatives 


Consider the global error E"+! = e"+! + RE", associated with the scheme 
(7). The transition operator 


A as oon 
R= (1 + Sey) (1 - Seu), 


is defined as follows: set z” = E, as the initial data with null boundary 
condition and zero source term f. After one time step of (7), let RE” be 


M 
the solution obtained. Using this, we have E?+! = S- R"-*e, 14. Thus, we 
k=0 
claim 
|R3|loo $C forall j=0,1,...,n. (8) 


Then it follows that sup ||E”*"||.. < C(As)?. Hence, the scheme (7) is 
nAs<s 


second-order accurate. It may be noted here that (8) is a stability condition. 
The details of this argument were given in [8, 8]. 


ot 
5th §) 
the local error with the method (7) satisfies 


<C forall (y, s)€9, for O0<i<3, then 


Lemma 2. If | 


lle] < C(As)°. (9) 


Proof. It can be shown using the same argument discussed in [3]. 


Theorem 2. The global error estimate E,, associated with (7) is given by 
|Z" loo < C(As)? for alln < S/As. 


Proof. See [8]. 


Theorem 3. The derivatives of z(y, s) satisfy the bounds 


oltm 


<Ce“"/?, (y, s)€D  foralll>0,m>0 with0<1+2m<8. 
Oy'ds™ 


Proof. Refer to [1] for the details. 
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3.4 Totally finite difference scheme 


The second-order approximation for the operator is given by 


n n n,n 
DtD-w? = — (3 mt “i1) 
Ye Ag + hppa Agta hy 


n-1 
In time, the backward difference scheme is D, 0; = , where wi = 


8 
w(y;, Sn). To solve (1), the Crank—Nicolson scheme for the time scale and the 
central difference scheme for the space are combined as 


2D, ZPt + LeZPtt = ptr ze Pt? _ ppze-? + peg fot — £.ZP (10) 
for 1<i<WM. 


Here, 
£.Z; = =eD SD, ZF + piZ?, f = (vis Sn), b} = b( yi, Sn), pi = p(yi)- 


After rearranging the terms in (10) and combining (7), the fully discrete 
scheme is 


1 

5 (sae +reZnth 4 eg =g!@, 1<i<M, 

aR = Gilead). Za = Oren), (11) 
Z,? = Ov(yi, —3;) for j=0,...,p and i=1<i<M, 


2 2 2 
ry =As(- =~), 1? = As( 2 + oP!) +1, rf =As(- = a ), 
h higths Aghia 


A 
gh = (—eptigpet —opap-r + eth + ap) + (1 Se.) ap 
for 0<i< M-1. 


The difference equation (11) at n +1 time level forms a tridiagonal system 
of M — 1 equations with the same number of unknowns. This system has 
properties like r; <0, r? > 0, ee <0 for 1<i2< M-—1. The Thomas 
algorithm [24] is used to solve the tridiagonal system. 


3.5 Error analysis 


Theorem 4. Let z be the analytical solution to (1) and let Z be the numer- 
ical solution (11). Then on S-mesh, the error of the scheme (11) satisfies the 
following estimate: 
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7 


ax |(z —Z) (yi, Sn)| < c((M—n My? + As”) for i=1,...,M-—1. 
On B-S-mesh, it satisfies 


max |(z — Z)(yi, sn)| <C(M-? + As?) for i=1,...,M—-1, 


a, 


where Z(yi, 8n) = ZP for (yi,8n) € D™. 


Proof. The proof is divided into various steps for each time level. On the 
first interval s € [0,7] (the time discretization n varies from 0 to p), the term 
f(y, s) —c(y, s)O4(s,s—7) in the right side of (1) is independent of «. Now 
since (y;, 81) € DM = beg x [0,-y), using the convergence results given in [3] 
and Theorem 2, we obtain on S-mesh 


max |(z— Z)(yi, 8n)| < c(M-?In? M+ (As)?) for O<i< WM. 


Following a similar procedure done in [831, 30] for the error bounds on B- 
S mesh, we consider the mesh generating function for B-S mesh ¢(€) = 
(2/B)Ex(s), where x(s) = —In (1 2(1 in)(s)). The mesh generating 
function satisfies ¢(€)’ < CM, and assume that « < M~!. The spatial 
mesh size h; on the layer region satisfies h; < CM —l and h; < Ce, for 
i = 0,1,...,(M/4) — 1 and for ¢ = (3M/4) + 1,...,M. Now using the 
bounds on the derivative of z given Theorem (3), we get on B-S-mesh, 


max |(z — Z)(yi, sn)| <C(M-? + As?) for i=1,...,M—-1, 


Now, the term z(y, s) depends on z(y, s — y), which is unknown for s > ¥. 
So, the above process is not applicable for s > y. To get the error over 
the interval [y, 2y], using the convergence results in [1] and Theorem 2, we 
obtain the desired bound. 


4 Post processing technique 


To enhance the order for the difference scheme (10), we use the Richardson 
extrapolation technique. First, we solve the discrete problems (11) on the 
fine mesh D2” = TI, x Tl.” with 2M and 2N mesh intervals in the spatial 
and time direction respectively, where i.” is the Shishkin-type mesh and 
is obtained by halving each mesh interval of TT” with a fixed transition 
parameter. Clearly, D” = {(yi, sn)} C D7” = {(9;,5n)}. Therefore, the 


corresponding S-mesh is i = {7,; € (0,1),0<i< 2M} by 
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‘ 
“ for 0<i< M/2, 
i(1—2 
7,= aS ME re le wee 
M 2 2 
i 
and the B-S-mesh is 
2 ; 
3( in (1-201 sh)ak) for 0<i<M/2-1, 
= YM/241 — YM/2-1\ /. M_...3M 
2 : M 
1 2( in (1-201 27) for ix 41<om, 


and 3, — 3,_1 = As/2 for 5; € TT”. Now, from Theorem 4, the error is 


(Z — z)(Yi, Sn) =c( M~'InM)/? + As?) +o((M~'In M)? + As”), 


( 
=c( 8) + CAs? + o((M- In M)? + As?) 


(12) 
for (yi, 8n) € D™. Let Z(G;, $n) be the solution to (11) on the domain D?™. 
From Theorem 4, we get 


FZ — 2)(U 8n) = C((2M)-*0Ge)? + (48°) + of (Mt In My? + (54)") 


(13 
for (Y;,5n) € D?™. Now, the elimination of o(M~?) term from (12) and (13) 
leads to the following approximation: 


Les 
2(yis 8n)—3 (42-2) (Yi, Sn) = o((M~'In M)?+As?), (Vis Sn) € D™. (14) 
Therefore, we use the extrapolation formula as 
Wa M 
Zestp(Yir 8n) = 3(42—Z) (yi, Sn), (Wis Sn) EDM. (15) 
Theorem 5. Let Zz) be the solution by extrapolation technique (15) and 
let z be the solution to problem (1). Also, assume that /é < M~1!. Then we 


have the following error bound on S-mesh 


2(Yi, Sn) Pa Zeatp (Yas Sn) 


<¢(M~‘In' M+ As‘) for tea. 


Proof. The term in right side of (1), f(y, s) — b(y, s)Ox(y, s — y) is inde- 
pendent of ¢ in the interval [0,7], where the time discretization parameter n 
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varies from 0 to p. One can find the complete analysis of the extrapolation 
technique in [28]. 


Since z(y, s) depends on z(y, s — y), which is unknown for s > ¥, the 
approach in [28] is not applicable in the interval [7,27] and s > 7. Therefore 
the following delay parabolic equation involving small parameter considered 
on domain D2 = II x (7, 27): 


(2 + Ley) 2(y, 8) = —b(y, 8)2(¥, 8-7) +H 8), (Ys 8) €Do, 


zy, s)= aly, s), (y, s)€ D1 =I x (0,4), (16) 
z(0, s) =@,(s), z(1, s)=0,(s), 5 € [7,29], 


where z,(y, 8) is a continuous solution in 1. Applying the scheme for (16) 
on D2 as given in (10), we have 


2D, Zt) — 62701 4 pi ZPtt = ppt gp eth _ pager (17) 
+h “tg aaa ~~ £-Z;, 
Zo = O1(Sn), ZN = On Sn) Sn € ly, 29], 
Z54 =Zr(yi, 8n), (yar 8n) € OM, 


where 67 = Dy Dy, f? = f(yis 8n), (Yi, 8n) € D4’ and Z;(-,-) is the numerical 
solution in D1". 


Decompose z in (16) on the domain Dz as z = w, + ws, where wy, is 
the smooth component and w, is the singular component. Again we write 
Wr = Wra + Wry, Satisfying 


3 + Pp(y)wroly, s) = —b(y, s)wroly, s— 7) +Fy, 8), (y, s) € Do, (18) 
wro(y, 8) — 2(y, 8), (y, s) ED, wro(0, s) = 2(0, s), se ly, 29], 


Owro YY, s) 
{ee 


and 


0 
(5, + Seu) weal, 8) = —b(y, 8) (¥,5— 7) + (wrodays (8) € D2, 
wri ly, s) — 0, (y, s) € D1, wr (0, 8) = 0, wri (1, 8) = 0, SE [y, 29]. 
(19) 
Therefore 


(2 + Ley )wrly, 8) = —b(y, s)wr(y, 8 = 7) ale iy, 8), (y, s) € Do, 


wrly, s) = 2(y, 8), (y, s) € D1, 
w,(0, s) = w,ro(0,8),  wr(1, s)=wro(l, 8), 8 € [7,24]. 
(20) 
Then, the singular component ws, satisfies 
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(2+ 2.,)wstu 8)=—b(y, s)rsly, 8). (Hr 8) € Da, 
ws(Y; s) 73 0, (y, s) € D1, w,(0, s) = 01 (s) — Wro (0, 8); 
w,(1, s) = 0,(s) — wro(1, s), SE [y, 27]. 


Write ws = Ws, + Ws,, where we, is the left boundary layer component on 0; 
and wes, is the right boundary layer component on 0,. Also, ws; and we, are 
satisfying the following PDEs: 


(ee soa s) =—b(y, s)wi(y, s—7), (y, 8) € Do, 


Os 
alt, s)= Wro (0, s) SE ee 24]; Wsi(Y, s) = 0, (y, s) € Du, 
as ae Le, ye ete s)= —b(y, s)wrly, s— pie (y, s) € Da, 
( 


ws, (1,8) = 8, — w-o(1, 8), 5 € [7,27], wsr(y, 5) =9, (y, 8) € Dir, 
(22) 
where 91; = (0,) x [0,7] and 91, = (1 — €,1) x [0,4]. Since, D2 C 9, the 
estimates given in Theorem 1, can be used for w, and w,. Decompose the 
numerical solution Z to (17) as Z = W,.+ W,, where W, is the smooth part 
and W, is the singular part. Thus 


n n+ n+i— n+3 
(Dy +D; +2.)W, Oe a We a arog eNO) 


We = Zo(Yis Sn), Ce Sn) € ait 
W,6 = w,(0, Sn), Wn = w,(1, Sin) Sn © 7; 2yI, 


and therefore, W, satisfies 


(Di +D; +2.)W, Wert? = —opW?t??, (yi,sn4n) DM, (24) 
Wet =0, (Yi, Sn) € DY", 
W.G = 01(sn) — wr(0, Sn), 
sv = 9Or(Sn) — wr(1, Sn), Sn € [y, 24). 


Now, we write W, = W., + Ws,, where W., is the boundary layer on 0)” and 
W,,. is the boundary layers on 9“. Hence W,, and W,, are defined by 


(pr D; £.)Warldi Sn41) — —bPWai(¥i, Spots), (Yi, Sn41) € Dy, 
Ws1(0, Sn) 7 ©; a wr (0, Sn) Cr Sn) e d1, Wily, Sn) —_ 0, (Yi, Sn) € ae 
(Di + Dy + £2)Werl¥is Sng) = —OPWorl¥is Saran) (ir Engg) © DM, 
W,,(1, Sn) = 0, ~~ wr (1, Sn)s (Yi, Sn) e o,, Ws, = 0, (Yi Si) € 3 
(25) 
where oY is the discretized domain of Dy; and D* is the discretized domain 
of D4}. Similarly OM is the discretized domain of s and Cea is the discretized 
domain of 0,. 
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4.1 Error bound for W, 


Lemma 3. The local truncation error in 93 associated to the smooth com- 
ponent satisfies 


(2D; + £) (wy —W,) (Yi, Sn) 
4 


— q—-1 2 2 iz : : 29 ac 1 
= c((M In M) + As ) + wit Yi-1) dy? (Yi, So) 
1 » Pw, —4 4 
pos Fee (Yi, 8n) + O(M™~* + As*). 


Proof. Using (20) and (23), we get 
(2D; + £2) (we —We)(gis Sn) 
= 88(2(Uis Spy) ~ ZolUir Sn—p-g)) + (2-205 )ewelis sn) 
+ Sworlti 8n—1) — +(Le,y — Le)Wr(Yis Sn—4). 


Now, by using the estimate in Theorem 4 and Taylor’s expansion, the desired 
result can be achieved. Refer to [8] for more details. 


Let the function Ea(y, s), for d = 1,2, satisfy the IBVPs (refer approach 
of Kellar [12]): 


(7) E Otw,(y, s) a 

— g = 
(5 + a) Ey (y, 8) 12 Oy! im D, (26) 
Ex(y, s) = 0, (y, s) € dp, 
E,(0, s)=0, Fi(1, s)=0, s€[0,S], 

a) 1 Bw,(y, s) 
(5. + Sou) Holy, )= pS ing, on 
Eo(y, s) — 0, (y, s) € Ts, 


E,(0, s)=0, £o(1, s)=0, s€[0,S]. 


Now Eq can be decomposed as Eq = na+Va, where ng and vq are the regular 
and singular layer parts of Eg. Now from Theorem 1, we have 


One can get these bounds using a similar procedure as done in [1]. Taking 
(y, 8) € De, (26) and (27) reduce to the following IBVPs: 


atm, 


si coat Cf en a a for 0<1+2m<8. 
Oy'ds™ 
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0 _& Otw,(y, 8). 
(= + L2,y) Erly, s)= 1 oy! in Do, 
Ex(y, s) = Ey, (y, s) € D1, 
E,(0, s) = 0, Ex(1, s) = 0, 8E ly, 29], 

) 1 Pw,(y, s) < 
C + £.,y) Ea(y, s)= ae in Do, 


Ey, s) = Foy, (y, s) € D1, 
E3(0, s) = 0, Fo(1, s) i 0, BE Ly, 2y], 


where E},(-,-) and £,(-,-) are the respective solutions in D,. Therefore, 
the components ng and Ug, d= 1,2, satisfy 


(2 +2n)n=Z 2 wo, 


Os | 12 Oy* 
0 1 Puw,(y, s) y 
Tr Ley) Me 9 Has in Do, 
naly, 8) = Eay, (y, s)€9D1 for d=1,2, 


na(0, s) = 0, na(1, s) = atl s), SE Be 24]; 


(2 +Ley)0a=0 ind, 


Valy, 8) =0, (y, s)€D, for d=1,2, (28) 
Va(0, 8) = 0, Va(1, s) 7 —na(1, 8), SE Ly, 27]. 


Lemma 4. The local truncation error in oY associated with w,. satisfies 


(wr — Wr) (yi; Sn) = c((uo} InM)? + As?) + (yi — yi—1)?m + mAs + O(M~* + As?). 
Proof. From Lemma 3 and (28), one can easily get 


(2D; a ) ((w, — W,)(yi; Sn) = C((M- In M)? + As) 


a((2 t Ley) 0; +29) n 


+As @ +4 i) -ODs4 2.) Ne 
+O(h* + As*). (29) 


Using the derivative bounds of ng and from the Taylor’s expansion, it follows 
that for d = 1, 2, 


a((2 5) (2D, 4 2))m va-((2 + Ley) — (2D; +29) no 
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< c(h! a As‘). (30) 
Therefore, from (29) and (30), we have 
(2D; ii 2, (wp —Wr)(Yiy Sn) < c((M—! In M)? + As? + (ht + As')), 


and, by applying barrier functions and the discrete maximum principle as 
done in [22], we obtain the following bound: 


(er — Weds Sn) 


<C (wr In M)? +As*) +h?m +n2As? + O(h* + As’). 


Lemma 5. The error associated with W,. after extrapolation satisfies 
(in = Wrevtp) We Sa) S c(M- + As‘) for (Yi, Sn) a. 


Proof. From Lemma 4 on the fine mesh 93%, we have 


2 


_ As?) h? As? 
(W, —wr) (yi, sn) = C((2M)~?(In2M)? ) bm +m 7 +0(M~4+4As*). (31) 


From the extrapolation formula (15), we can write 
—— 
(wy — W,-eatp) (Yis Sn) = Wr(Yi, $n) — (S4W, —W,)(yi, 8n)) 
1 — 
= —5 (401; — w,) + (We — wr) is $n): 


By using Lemma 4 and (31), we obtain 


-; (40, =) +(We= wr)) (yi, Sn) =; ( = c(i In2M)? + As?) 


= hem = mAs” 
*: c(i In M)? + As?) 


+ h?m + mae) + O(M-* + As*) 


=O0(M-* + As‘), 


which is our desired bound. 


Now, we define the function Fy = Fa + Far, d = 1,2, by the following 
IBVPs: 
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(2 oz Ley) Fu = _ pot O wel. 8) (y, s) in D; = (0, g) as (0, S], 


Os 3 Oy* 
) Ows(y, ; 
(Z+2ey)Fir= EAM) Gy, 5) in D. = (1-§1) x 0.5}, 
Fuly, 8) =0, (y, s) € [0,€] x (7,9), 
Fir(y, s)=0, (y, sje [mee 1) x (e a) 
F,(0, 8s) = Fii(€, s)=0, SE (0, S], 
Fy,(1— €, s) = F,,(1, s)=0, SE (0, S], 
(32) 
7) OPwesly, 8 
(= a Ley) Fa ja a (y, s) Mn Di = (0, €) x (0, S], 
2) Pus(y, s 
(= 4: Ley) Far = a (y, s) in D, = (1-€,1) x (0, S], 
Frly, 8) =0, (y, 8) € [0,€] x (7,0), 
F2,(y, s) =0, (y, s)€ [LSS 1] x ic 7,0), 
F(0, s)= Fiil€, s)= =0, SE (0,S], 
Fy,(1— €, s)= F,,(, s)=0, SE (0, S], 
(33) 
The solution Fy, d = 1,2, to (32) and (33) satisfies the following bounds: 
gltm 
Saat] <Ce?(expl—yV/ B/E) + ex(—— VBP), us 5) 2. 


In this context, by considering s € (y, 27), therefore (32) and (33) reduces to 


) Ows(y, : 
(= a Ley) Fu ras See (y, s) in D) = (0, €) x (7,27), 
) Ows(y, 
(= 4+ ey) Fir = EO al (y, s) im Dy, = (1—€,1) x (y, 27); 
Fu(y, 8) = Fin ly, 8), (y, s) € (0, €] x (0,7), 
Fir(y, 8) = Firy(y, §), (y, s) €[1—€, 1] x (0,7), 
F(0, s) = Fu(€,t) =0, SE (7,27), 
F,.(1 — €,t) = Fi,(1, s) =0, SE (y, 24), 


(34) 
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Bau. (y, ; 
(2 ne Lex) Fa = A OO) (y, 8) in D1 = (0,€) x (7,27), 


Os 19 0s? 

O 1 Pus(y, : 
G oy Ley) For ~ 19 UES (y, 8) in Dp = (1—€,1) x (7,2); 
Fu(y, 8) = Fay(y, a) (y, s) € (0, €] x (0,7), 
Fo,(y, 8) cae For (y, 8), (y, s) € [1 as g, 1 x (0,7), 
F2(0, 8) = Fy (€, s)= ’ s € (7,29), 
Fo, (1 — €, s) = F,,(1, s) =0, SE (7, 27), 

(35) 


where Fyi(-,-), k = 1,2, are the analytic solutions in [0, €] x (0, y) and Fy,(-,-), 
k = 1,2, are the analytic solutions in [1 — €,1] x (0,7). 


Lemma 6. For (y;, 5,) D2’, the error associated with W, satisfies 
(W.—ws) (Yi, 8n) = (M7" In M)? Fy (y;, 8n)+As? Fo(y:, sn)+O(M~4 In4 M+As*). 


Proof. Write W, = Ws, +W;,, where W,, and W,,. are defined in (25). Now 
W,—ws = (Ws,—Ws1) +(Ws, —Ws,), where the error W,;— ws, is related to a 
layer on D; and W,,.— ws, on D,. These errors can be estimated separately. 
First, we estimate W,, — ws,, from (22) and difference equation (25) and 
follow the similar procedure of Lemma 4. 


Next, we can estimate the error W;,. — Ws,. 


Lemma 7. For (yi, Sn) ®,™ , the error associated to W, after extrapolation 
satisfies 
(ws — Weeatp) (Yas Sn) 


< c(M-*Int M+ As‘). 
Proof. From Lemma 6, we get 
(W,—ws)= (M-!InM)?F\(y, sn) + As? Fo(y%, Sn) 


+0(M~4 In* N + As‘) fOr is 85 Oo, (36) 


and 


— A 2 
(ws — We) = ((2M)~" m2M)?Fi(ys, 8) + Fay, $n) 


+0 (4 Int M+ As‘) (37) 


for (yi, Sn) € D3™. Eliminating the terms O(M~?) and As? from (36) and 
(37), the required result is achieved. 


Theorem 6 (Error after extrapolation in D2‘). Let Zertp be the extrapo- 
lated solution (by technique (15)) for solving (11) on D,™ and D3”. Let z 
be the solution to (1). Assume that ¢ < M~?. Then we have the following 
error bound associated with Zeztp : 
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2(Yi, Sn) — Lexty Yes Sn) 


<¢(M~‘In' M+ As‘) for (“heSS A: 


Proof. We have 


Z(yi, Sn) — Zextp(Yis Sn) < Wr(Yi, Sn) = Wy eatp(Yis Sn) 


i Ws(Yis 8n) — Weeatp(Yis Sn) 


for all (yi, Sn) € D3’. Combining Lemmas 4 and 7, we can get the required 
result. 


Remark 1. The error bound on B-S-mesh in D3” is 


2(Yi, a eee Sr mE Sn) <c(M~‘ + As‘) for 1 <i< M-1. 
Remark 2.(Error after extrapolation in D”) Let Zextp be the extrapo- 

lated solution (by technique (15)) for solving (11) on D,™ and D3”. Let z 

be the solution to (1). Then we have the following error bound on S-mesh: 


Z(Yi, Sn) — Zeatp(Yi; Sn)| < c(M~ Int M+ As”) for 1<i<M-1, 


similarly, on B-S-mesh 


2(Yi, Sn) 2 Leseg Uys Sn) 


<c(M~* + As‘) for 1<i<M-—1l. 


5 Numerical experiments 


The proposed scheme (10) is tested on two test problems in this section. In 
this section, all the numerical results are obtained using MATLAB software 
in 64 GB RAM workstation. 


Example 1. Consider 


es EZyy + 0.52 = —2e~tz(y, a 1) Tr f(y, 8), (x, t) € (0, 1) x (0, 2], 
Bi Bie PUNE ee OTE 8) E1051) oT, 0) 
2(0, 8s) =e? +e St VV 2(1, 8) =e St Ve +e-*, 5 € (0, QI. 

(38) 


The exact solution for Example 1 is z(y, s) = e~8+¥/V@ + e~stC-w)/véE, 
The maximum pointwise error before and after extrapolation given by 


Ble = mae. Us ta) Ze an) 


(Yir Sn)EDM 
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and 


M,As M,As 
E! =Z 


= je gia 0 Sn) extp (Yi Sn)|. 


. . EM.As 
The corresponding order of convergence is P-4* = log, pr ae7s |- Here, 
€ 
M,As 


2(yi, Sn) is the exact solution and Z-4%(y;, sn) and Zerip (Yi, Sn) are the 
numerical solutions before and after extrapolation, respectively. 


Example 2. Consider the test problem: 
1 2 
Bg = Egy F ( wy z=s?—2z(y,s—1), (y, s) € (0, 1) x (0, 2], 


z(y, s) = 0, (y, s) € (0, 1] x [=1; 0], 
2(0,s) =0,2(1,s) =0, s € [0, 2). 


(39) 


Since the exact solution to (2) is not known, we use the idea of double 
mesh principle to obtain the pointwise errors and to verify the ¢-uniform 
convergence. Define Z(y;, sn) as the numerical solution obtained on 9?” = 
Teé x 12% with 2M mesh intervals in space and 2N mesh intervals in the 
s-direction, where ie is the Shishkin-type mesh as defined II)” with the 
fixed transition parameter. For each ¢, we calculate the maximum pointwise 


error before and after extrapolation by E)4* = ae |Z A5(y;, Sn)— 
Yir 8n)ED 

7M,As GM,As _ M,As 7M,As 

Z (yi, Sn)| and E: ee max [erp (Yi, Sn) _ Desin (Yi; Sn) |, 


(yi ; Sn)ED MM. 
respectively. The corresponding order of convergence is obtained by P”4* = 


EM As >M,A . : ‘ . 
logs grinasys |- Here, Zerip (Yi, Sn) is the extrapolation solution obtained 
€ 


by the double mesh principle. 


n 0.5 me 


No 


Figure 1: Surface plots of the computed solution on S-mesh for Example 1. 
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e=10? 


t-0.5 


Figure 2: 
e=10% 
x10 
4 
nN 0.5 
0 
2 
‘ 0.5 
3 0 0 y 


(a) Before extrapolation. 


e=10-° 


Solution plots of each time on S-mesh for Example 1. 


(b) After extrapolation. 


Figure 3: Error graphs of the computed solution on B-S mesh for Example 1. 
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(a) Before extrapolation. 


Surface plots of the computed solution on S-mesh for Example 2. 
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Maxx Error 
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(b) After extrapolation. 


Figure 9: Log-log plots on B-S mesh for Example 2. 
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e=10° 


0 2 
19 SK = 107 10 <= 1072 
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(a) Before extrapolation. 


(b) After extrapolation. 


Figure 6: Log-log plots on S-mesh for Example 1. 


1 
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$ 

103 = 
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(a) Before extrapolation. 


(b) After extrapolation. 


Figure 7: Log-log plots on B-S mesh for Example 1. 
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(a) Before extrapolation. 
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(b) After extrapolation. 


Figure 8: Log-log plots on S-mesh for Example 2. 


Table 1: £2°4° and P“4* generated on S-mesh using the proposed method for 


Example 1. 
M/As_ | Extrapolation le-4 le-6 
Singular layer | Regular layer | Singular layer | Regular layer 
32/10 before 1.15e-2 6.42e-4 1.15e-2 6.42e-4 
1.47 2.13 1.47 2.13 
after 3.09e-4 1.52e-6 3.09e-4 1.52e-6 
2.72 3.56 2.72 3.56 
64/40 before 4.15e-3 1.46e-4 4.15e-3 1.46e-4 
56 2.25 1.56 2.25 
after 1.33¢-6 1.10e-7 1.33e-6 1.10e-7 
3.05 3.78 3.05 3.78 
128/160 before 1.40e-3 3.07e-5 1.40e-3 3.07e-5 
65 2.38 1.65 2.38 
after 5.61e-6 9.33e-9 5.61e-6 9.33e-9 
3.22 3.89 3.22 3.89 
256/640 before 4.46¢e-4 5.90e-6 4.46¢e-4 5.90e-6 
1.68 2.48 1.68 2.48 
after 6.01e-7 6.28e-10 6.00e-7 6.28e-10 
3.34 4.27 3.34 4.27 
512/2560 before 1.38¢-4 1.05e-6 1.38¢-4 1.05e-6 
after 5.91e-8 3.25e-11 5.91e-8 3.25e-11 
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Table 2: EM As and pMAs generated on B-S-mesh using the proposed method for 
Example 1. 
E Extrapolation Number of intervals 
32/10 64/40 128/160 256/640 | 512/2560 
le-—6 before 4.9562e-3 | 1.3618e-3 | 3.4560e-4 | 8.6895e-5 | 2.1776e-5 
1.8637 1.9783 1.9918 1.9965 
after 5.3442e-4 | 4.4861e-5 | 3.0539e-6 | 1.9448e-7 | 1.2234e-8 
3.5744 3.8767 3.9730 3.9906 
le—8 before 4.9562e-3 | 1.3618e-3 | 3.4560e-4 | 8.6895e-5 | 2.1776¢-5 
1.8637 1.9783 1.9918 1.9965 
after 5.3442e-4 | 4.4861e-5 | 3.0539e-6 | 1.9448e-7 | 1.2234e-8 
3.5744 3.8767 3.9730 3.9906 
le — 12 before 4.9562e-3 | 1.3618e-3 | 3.4560e-4 | 8.6895e-5 | 2.1776e-5 
1.8637 1.9783 1.9918 1.9965 
after 5.3442e-4 | 4.4861e-5 | 3.0539e-6 | 1.9448e-7 | 1.2234e-8 
3.5744 3.8767 3.9730 3.9906 


Table 3: 624° and P¥4s generated on S-mesh using the proposed method for 


Example 2. 
N/As Extrapolation le-4 le-6 le-8 
Inner layer | Outer layer | Inner layer | Outer layer | Inner layer | Outer layer 
32/10 efore 9.021e-3 2.800¢-3 9.021e-3 3.100e-3 9.021e-3 3.100e-3 
633 2.001 418 2.015 1.418 2.015 
after 1.850e-5 6.976¢e-6 3.344e-4 7.103e-6 3.344e-4 7.103e-6 
3.788 3.982 2.704 3.970 2.704 3.970 
64/40 efore 2.907e-3 7.068¢-4 3.374e-3 7.669¢e-4 3.374e-3 7.669e-4 
911 2.001 528 2.000 1.528 2.000 
after 1.339e-6 4.414e-7 5.132e-5 4.530e-7 5.132e-5 4.530¢e-7 
3.878 4.06 3.037 3.925 3.037 3.925 
128/160 efore 7.729¢-4 -765e-4 1.169e-3 -917e-4 1.169e-3 -917e-4 
418 2.002 633 2.000 1.6337 2.000 
after 9.108e-8 2.631e-8 6.249e-6 2.981e-8 6.249e-6 2.981e-8 
3.949 4.024 3.489 4.011 3.489 4.011 
256/640 efore 1.959e-4 4.404e-5 3.979e-4 4.793¢e-5 3.979e-4 4.793¢e-5 
.980 2.010 655 2.000 1.655 2.000 
after 5.897e-9 617e-9 5.565e-7 .848¢e-9 5.565e-7 .848e-9 
3.974 4.162 3.920 4.291 3.920 4.291 
512/2560 efore 4.909e-5 -093e-5 1.263e-4 -198e-5 1.263e-4 -198e-5 
after 3.752e-10 1.778e-11 3.675e-8 9.440e-11 3.675e-8 9.440e-11 
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Table 4: pias and pias generated on B-S-mesh using the proposed method for 
Example 2. 
E Extrapolation Number of intervals 
32/10 64/40 128/160 256/640 | 512/2560 
le-—6 before 1.3012¢e-2 | 3.5208e-3 | 9.2402e-4 | 2.3350e-4 | 5.8667e-5 
1.8859 1.9299 1.9845 1.9928 
after 5.4243e-4 | 4.6020e-5 | 3.1547e-6 | 2.0324e-7 | 1.2801e-8 
3.5591 3.8667 3.9562 3.9889 
le—8 before 1.3012e-2 | 3.5208e-3 | 9.2402e-4 | 2.3350e-4 | 5.8667e-5 
1.8859 1.9299 1.9845 1.9928 
after 5.4243e-4 | 4.6020e-5 | 3.1547e-6 | 2.0324e-7 | 1.2801e-8 
3.5591 3.8667 3.9562 3.9889 
le — 12 before 1.3012e-2 | 3.5208e-3 | 9.2402e-4 | 2.3350e-4 | 5.8667e-5 
1.8859 1.9299 1.9845 1.9928 
after 5.4243e-4 | 4.6020e-5 | 3.1547e-6 | 2.0324e-7 | 1.2801e-8 
3.5591 3.8667 3.9562 3.9889 
The surface plots are plotted for ¢ = 107? and ¢ = 10~® in Figure 1 for 


Example 1 and Figure 4 for Example 2 on S-mesh, respectively. Figures 2 and 
5 display the solution for different values ¢ with respect to time for Example 
1 and Example 2, respectively. From these figures, one can visualize that the 
boundary layers appear at y = 0 and y = 1. The error plots are given in 
Figure 3 for Example 1 on the B-S mesh before and after extrapolation. These 
graphs show that the error is less after extrapolation compared to before 
extrapolation. To see the numerical convergence rate graphically, E4° are 
plotted in the log-log scale before extrapolation in Figures 6(a), 7(a),8(a), and 
9(a) and after extrapolation in Figures 6(b), 7(b),8(b), and 9(b). Moreover, 
Bes and P ‘AS by the proposed scheme for Example 1 are presented in 
Tables 1 and 2 on S-mesh and B-S-mesh, respectively. Similar results for 
Example 2 are shown in Tables 3 and 4. From these tables, one can conclude 
that B-S-meshes are more accurate, and the rate of convergence is more 
compared to the S-mesh. One can notice that the numerical results are well 
by our theoretical findings. 
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